############################################################
## APPENDIX REPRODUCTION SCRIPT (A3–A6, A8–A11)
## Inputs:
##   - merged_experiment_data.sav
##   - merged_USwaves_for_FigA9.sav
## Outputs (standardized):
##   A3:  A3_Table_sample_composition.csv
##   A4:  A4_Table_region_effects_combined.csv
##   A5:  A5_1_Table_H1_pooled.csv
##        A5_2_Table_H1_by_country.csv
##   A6:  A6_1_Table_H2_H3_H4.csv
##   A8:  A8_1_Table_all_treatment_effects.csv
##        A8_2_Table_context_pooled.csv
##        A8_3_Table_context_by_country.csv
##        A8_Fig1_contextual_factors_coefplot.png
##   A9:  A9_Table_USwaves_full.csv
##        A9_Table_USwaves_ACpass.csv
##        A9_Fig_data_party3_by_ally_and_wave.csv
##        A9_Fig_party3_by_ally_and_wave.png
##  A10:  A10_Table_power_curve_grid_ANALYTICAL.csv
##        A10_Table_MDE_80pct_grid_ANALYTICAL.csv
##        A10_Fig_power_curves_ANALYTICAL.png
##  A11:  A11_Table_meta_inputs_ally_country_estimates.csv
##        A11_Table_meta_random_effects_summary.csv
##        A11_Table_meta_overall_OR.csv
##        A11_Fig_meta_ally_forestplot.png
############################################################

suppressPackageStartupMessages({
  library(haven)
  library(dplyr)
  library(tidyr)
  library(MASS)      # polr
  library(tibble)
  library(ggplot2)
  library(broom)     # meta block
  library(metafor)   # meta block
})

# ============================================================
# 0) Global setup
# ============================================================

countries <- c("Brazil", "China", "India", "Russia", "UK", "USA")

base_family <- "Garamond"
if (.Platform$OS.type == "windows") {
  suppressWarnings(try(windowsFonts(Garamond = windowsFont("Garamond")), silent = TRUE))
  if (!"Garamond" %in% names(windowsFonts())) base_family <- "serif"
}

theme_pub <- ggplot2::theme_minimal(base_size = 15, base_family = base_family) +
  ggplot2::theme(
    text = ggplot2::element_text(family = base_family),
    panel.grid.minor = ggplot2::element_blank()
  )

forest_style <- function() {
  theme_pub +
    ggplot2::theme(
      panel.grid.major.y = ggplot2::element_blank(),
      panel.grid.minor   = ggplot2::element_blank(),
      axis.title.y       = ggplot2::element_blank()
    )
}

data_all <- haven::read_sav("merged_experiment_data.sav") %>%
  dplyr::mutate(
    dataset = as.character(dataset),
    strike5 = factor(as.character(strike5),
                     levels = as.character(0:4),
                     ordered = TRUE)
  )

data_clean <- data_all %>% dplyr::filter(dataset != "Original")

# ============================================================
# 1) Helpers
# ============================================================

tidy_polr_H <- function(fit) {
  N_fit <- stats::nobs(fit)
  cs <- coef(summary(fit))
  td <- as.data.frame(cs)
  td$term <- rownames(td)
  rownames(td) <- NULL
  
  td %>%
    dplyr::rename(
      estimate  = Value,
      std.error = `Std. Error`,
      statistic = `t value`
    ) %>%
    dplyr::mutate(
      conf.low  = estimate - 1.96 * std.error,
      conf.high = estimate + 1.96 * std.error,
      p.value   = 2 * (1 - stats::pnorm(abs(statistic))),
      N         = N_fit
    ) %>%
    dplyr::filter(!grepl("\\|", term))
}

label_predictor <- function(term) {
  dplyr::case_when(
    term == "ally" ~ "Alliance (vs no alliance)",
    term == "male" ~ "Male (vs female)",
    term == "age"  ~ "Age (in decades)",
    term == "ed4"  ~ "College+ (vs less than college)",
    term == "dem"  ~ "Democracy (vs non-democracy)",
    term == "ints" ~ "High stakes (vs low stakes)",
    term == "cost" ~ "Low cost (vs high cost)",
    grepl("^region_exp", term) ~ {
      lvl <- gsub("^region_exp", "", term)
      paste0("Region ", lvl, " (vs 1)")
    },
    grepl("^dataset", term) ~ {
      lvl <- gsub("^dataset", "", term)
      paste0(lvl, " (vs USA)")
    },
    TRUE ~ term
  )
}

# ============================================================
# A3 — Table A3 (Sample composition)
# ============================================================

data_A3 <- data_all %>%
  dplyr::filter(dataset %in% countries) %>%
  dplyr::mutate(
    age_years = as.numeric(age) * 10,
    age_group = cut(
      age_years,
      breaks = c(18, 25, 30, 40, 50, 60, 70, Inf),
      right  = FALSE,
      labels = c("18–24", "25–29", "30–39", "40–49",
                 "50–59", "60–69", "70+")
    )
  )

N_long <- data_A3 %>%
  dplyr::group_by(dataset) %>%
  dplyr::summarise(value = dplyr::n(), .groups = "drop") %>%
  dplyr::mutate(Variable = "N", Category = "")

gender_long <- data_A3 %>%
  dplyr::filter(!is.na(male)) %>%
  dplyr::mutate(
    Variable = "Gender",
    Category = ifelse(as.numeric(male) == 1, "Male", "Female")
  ) %>%
  dplyr::count(dataset, Variable, Category, name = "n") %>%
  dplyr::group_by(dataset, Variable) %>%
  dplyr::mutate(value = 100 * n / sum(n)) %>%
  dplyr::ungroup() %>%
  dplyr::select(dataset, Variable, Category, value)

educ_long <- data_A3 %>%
  dplyr::filter(!is.na(ed4)) %>%
  dplyr::mutate(
    Variable = "Education",
    Category = ifelse(as.numeric(ed4) == 1, "College+", "Less than college")
  ) %>%
  dplyr::count(dataset, Variable, Category, name = "n") %>%
  dplyr::group_by(dataset, Variable) %>%
  dplyr::mutate(value = 100 * n / sum(n)) %>%
  dplyr::ungroup() %>%
  dplyr::select(dataset, Variable, Category, value)

age_long <- data_A3 %>%
  dplyr::filter(!is.na(age_group)) %>%
  dplyr::mutate(
    Variable = "Age",
    Category = as.character(age_group)
  ) %>%
  dplyr::count(dataset, Variable, Category, name = "n") %>%
  dplyr::group_by(dataset, Variable) %>%
  dplyr::mutate(value = 100 * n / sum(n)) %>%
  dplyr::ungroup() %>%
  dplyr::select(dataset, Variable, Category, value)

A3_Table <- dplyr::bind_rows(
  N_long %>% dplyr::select(dataset, Variable, Category, value),
  gender_long, age_long, educ_long
) %>%
  dplyr::mutate(
    value    = ifelse(Variable == "N", value, round(value, 1)),
    Variable = factor(Variable, levels = c("N", "Gender", "Age", "Education"))
  ) %>%
  dplyr::arrange(Variable, Category) %>%
  dplyr::select(Variable, Category, dataset, value) %>%
  tidyr::pivot_wider(names_from = dataset, values_from = value) %>%
  dplyr::arrange(Variable, Category)

utils::write.csv(A3_Table, "A3_Table_sample_composition.csv", row.names = FALSE)

# ============================================================
# A4 — Table A4 (Region robustness; identical logic)
# ============================================================

region_results <- list()

for (nm in countries) {
  df_nm <- data_all %>%
    dplyr::filter(dataset == nm) %>%
    dplyr::filter(!is.na(region_exp)) %>%
    dplyr::mutate(region_exp = factor(as.character(region_exp)))
  
  if ("1" %in% levels(df_nm$region_exp)) {
    df_nm <- df_nm %>% dplyr::mutate(region_exp = stats::relevel(region_exp, ref = "1"))
  }
  
  if (nrow(df_nm) < 10L || length(unique(df_nm$region_exp)) < 2L) next
  
  fit_region <- MASS::polr(
    strike5 ~ ally + dem + ints + cost + region_exp + male + age + ed4,
    data = df_nm, Hess = TRUE
  )
  
  region_results[[nm]] <- tidy_polr_H(fit_region) %>%
    dplyr::mutate(dataset = nm, Predictor = label_predictor(term), .before = 1) %>%
    dplyr::select(dataset, Predictor, estimate, std.error, p.value, N, term)
}

A4_Table <- dplyr::bind_rows(region_results) %>%
  dplyr::mutate(dataset = factor(dataset, levels = countries)) %>%
  dplyr::arrange(dataset, Predictor)

utils::write.csv(A4_Table, "A4_Table_region_effects_combined.csv", row.names = FALSE)

# ============================================================
# A5 — Tables A5.1 (H1 pooled) and A5.2 (H1 by country)
# ============================================================

data_H1_pooled <- data_all %>%
  dplyr::filter(dataset != "Original") %>%
  dplyr::mutate(dataset = stats::relevel(factor(dataset), ref = "USA"))

fit_H1_pooled <- MASS::polr(
  strike5 ~ ally + male + age + ed4 + dataset,
  data = data_H1_pooled, Hess = TRUE
)

A5_1_Table <- tidy_polr_H(fit_H1_pooled) %>%
  dplyr::mutate(dataset = "Pooled", Predictor = label_predictor(term), .before = 1) %>%
  dplyr::select(dataset, Predictor, estimate, std.error, p.value, N, term)

utils::write.csv(A5_1_Table, "A5_1_Table_H1_pooled.csv", row.names = FALSE)

country_tabs <- list()
for (cty in countries) {
  df_cty <- data_H1_pooled %>% dplyr::filter(as.character(dataset) == cty)
  
  fit_cty <- MASS::polr(
    strike5 ~ ally + male + age + ed4,
    data = df_cty, Hess = TRUE
  )
  
  country_tabs[[cty]] <- tidy_polr_H(fit_cty) %>%
    dplyr::mutate(dataset = cty, Predictor = label_predictor(term), .before = 1) %>%
    dplyr::select(dataset, Predictor, estimate, std.error, p.value, N, term)
}

A5_2_Table <- dplyr::bind_rows(country_tabs) %>%
  dplyr::mutate(dataset = factor(dataset, levels = countries)) %>%
  dplyr::arrange(dataset, Predictor)

utils::write.csv(A5_2_Table, "A5_2_Table_H1_by_country.csv", row.names = FALSE)

# ============================================================
# A6 — Table A6.1 (H2–H4 interaction tests)
# ============================================================

data_H2 <- data_clean %>%
  dplyr::filter(dataset %in% c("USA","UK","Russia","India","Brazil","China")) %>%
  dplyr::mutate(
    NATO = ifelse(dataset %in% c("USA","UK"), "NATO", "Non-NATO"),
    NATO = factor(NATO, levels = c("Non-NATO","NATO"))
  )

fit_H2 <- MASS::polr(strike5 ~ ally * NATO + male + age + ed4, data = data_H2, Hess = TRUE)

data_H3 <- data_clean %>%
  dplyr::filter(dataset %in% c("Russia","China","India","Brazil")) %>%
  dplyr::mutate(
    Rus_vs_CIB = ifelse(dataset == "Russia", "Russia", "China/India/Brazil"),
    Rus_vs_CIB = factor(Rus_vs_CIB, levels = c("China/India/Brazil","Russia"))
  )

fit_H3 <- MASS::polr(strike5 ~ ally * Rus_vs_CIB + male + age + ed4, data = data_H3, Hess = TRUE)

data_H4 <- data_clean %>%
  dplyr::filter(dataset %in% c("USA","UK")) %>%
  dplyr::mutate(US_vs_UK = factor(dataset, levels = c("USA","UK")))

fit_H4 <- MASS::polr(strike5 ~ ally * US_vs_UK + male + age + ed4, data = data_H4, Hess = TRUE)

get_stats <- function(fit) tidy_polr_H(fit) %>% dplyr::select(term, estimate, std.error, p.value)

tab_H2 <- get_stats(fit_H2); N2 <- stats::nobs(fit_H2)
tab_H3 <- get_stats(fit_H3); N3 <- stats::nobs(fit_H3)
tab_H4 <- get_stats(fit_H4); N4 <- stats::nobs(fit_H4)

rows_A61 <- tibble::tibble(
  row_order = 1:11,
  Predictor = c(
    "NATO (vs non-NATO)",
    "Alliance × NATO",
    "Russia (vs China / India / Brazil)",
    "Alliance × Russia",
    "UK (vs USA)",
    "Alliance × UK",
    "Alliance (vs no alliance)",
    "Male (vs female)",
    "Age (in decades)",
    "College+ (vs less than college)",
    "N"
  ),
  H2_term = c("NATONATO", "ally:NATONATO", NA, NA, NA, NA, "ally", "male", "age", "ed4", NA),
  H3_term = c(NA, NA, "Rus_vs_CIBRussia", "ally:Rus_vs_CIBRussia", NA, NA, "ally", "male", "age", "ed4", NA),
  H4_term = c(NA, NA, NA, NA, "US_vs_UKUK", "ally:US_vs_UKUK", "ally", "male", "age", "ed4", NA)
)

pull_model_cols <- function(rows, tab, term_col, prefix, Nfit) {
  rows %>%
    dplyr::transmute(row_order, Predictor, term_name = .data[[term_col]]) %>%
    dplyr::left_join(
      tab %>% dplyr::mutate(term_name = term) %>% dplyr::select(term_name, estimate, std.error, p.value),
      by = "term_name"
    ) %>%
    dplyr::select(row_order, Predictor, estimate, std.error, p.value) %>%
    dplyr::rename_with(~ paste0(prefix, "_", .x), c(estimate, std.error, p.value)) %>%
    dplyr::mutate(!!paste0(prefix, "_N") := ifelse(Predictor == "N", Nfit, NA_real_)) %>%
    dplyr::select(-dplyr::all_of("Predictor"))
}

A61_H2 <- pull_model_cols(rows_A61, tab_H2, "H2_term", "H2", N2)
A61_H3 <- pull_model_cols(rows_A61, tab_H3, "H3_term", "H3", N3)
A61_H4 <- pull_model_cols(rows_A61, tab_H4, "H4_term", "H4", N4)

A6_1_Table <- rows_A61 %>%
  dplyr::select(row_order, Predictor) %>%
  dplyr::left_join(A61_H2, by = "row_order") %>%
  dplyr::left_join(A61_H3, by = "row_order") %>%
  dplyr::left_join(A61_H4, by = "row_order") %>%
  dplyr::arrange(row_order) %>%
  dplyr::select(-row_order)

utils::write.csv(A6_1_Table, "A6_1_Table_H2_H3_H4.csv", row.names = FALSE)

# ============================================================
# A8 — H6 contextual factors
# ============================================================

data_no_original <- data_all %>% dplyr::filter(dataset != "Original")

data_A81_pooled <- data_no_original %>%
  dplyr::mutate(dataset = stats::relevel(factor(dataset), ref = "USA"))

fit_A81_pooled <- MASS::polr(
  strike5 ~ ally + dem + ints + cost + male + age + ed4 + dataset,
  data = data_A81_pooled, Hess = TRUE
)

tab_A81_pooled <- tidy_polr_H(fit_A81_pooled) %>% dplyr::select(term, estimate, std.error, p.value)
N_pooled <- stats::nobs(fit_A81_pooled)

tabs_A81_country <- list()
Ns_country <- list()

for (cty in countries) {
  df_cty <- data_no_original %>% dplyr::filter(dataset == cty)
  
  fit_cty <- MASS::polr(
    strike5 ~ ally + dem + ints + cost + male + age + ed4,
    data = df_cty, Hess = TRUE
  )
  
  tabs_A81_country[[cty]] <- tidy_polr_H(fit_cty) %>% dplyr::select(term, estimate, std.error, p.value)
  Ns_country[[cty]] <- stats::nobs(fit_cty)
}

rows_A81 <- tibble::tibble(
  row_order = 1:13,
  Effect = c(
    "Alliance (vs no alliance)",
    "Democracy (vs non-democracy)",
    "High stakes (vs low stakes)",
    "Low cost (vs high cost)",
    "Male (vs female)",
    "Age (in decades)",
    "College+ (vs less than college)",
    "Brazil (vs USA)",
    "China (vs USA)",
    "India (vs USA)",
    "Russia (vs USA)",
    "UK (vs USA)",
    "N"
  ),
  term = c(
    "ally", "dem", "ints", "cost", "male", "age", "ed4",
    "datasetBrazil", "datasetChina", "datasetIndia", "datasetRussia", "datasetUK",
    NA_character_
  )
)

make_cols <- function(rows, tab, prefix, Nfit) {
  rows %>%
    dplyr::left_join(tab, by = "term") %>%
    dplyr::select(row_order, estimate, std.error, p.value) %>%
    dplyr::rename_with(~ paste0(prefix, "_", .x), c(estimate, std.error, p.value)) %>%
    dplyr::mutate(!!paste0(prefix, "_N") := ifelse(rows$Effect == "N", Nfit, NA_real_))
}

A81_pooled_cols <- make_cols(rows_A81, tab_A81_pooled, "Pooled", N_pooled)
A81_country_cols <- lapply(countries, function(cty) make_cols(rows_A81, tabs_A81_country[[cty]], cty, Ns_country[[cty]]))

A8_1_Table <- rows_A81 %>%
  dplyr::select(row_order, Effect) %>%
  dplyr::left_join(A81_pooled_cols, by = "row_order")

for (j in seq_along(countries)) {
  A8_1_Table <- A8_1_Table %>% dplyr::left_join(A81_country_cols[[j]], by = "row_order")
}

A8_1_Table <- A8_1_Table %>% dplyr::arrange(row_order) %>% dplyr::select(-row_order)
utils::write.csv(A8_1_Table, "A8_1_Table_all_treatment_effects.csv", row.names = FALSE)

fig_A81_df <- tidy_polr_H(fit_A81_pooled) %>%
  dplyr::filter(term %in% c("ally", "dem", "ints", "cost")) %>%
  dplyr::mutate(
    Effect = dplyr::case_when(
      term == "ally" ~ "Alliance (vs no alliance)",
      term == "dem"  ~ "Democracy (vs non-democracy)",
      term == "ints" ~ "High stakes (vs low stakes)",
      term == "cost" ~ "Low cost (vs high cost)",
      TRUE ~ term
    ),
    Effect = factor(Effect, levels = c(
      "Low cost (vs high cost)",
      "High stakes (vs low stakes)",
      "Democracy (vs non-democracy)",
      "Alliance (vs no alliance)"
    ))
  )

A8_Fig1 <- ggplot(fig_A81_df, aes(y = Effect, x = estimate)) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), height = 0.15) +
  geom_point(size = 2) +
  labs(
    title = "Appendix 8 Figure 1: Contextual factors (pooled)",
    x = "Estimate (log-odds, ordinal logit)",
    y = NULL
  ) +
  theme_pub

ggsave("A8_Fig1_contextual_factors_coefplot.png", A8_Fig1, width = 7.5, height = 4.8, dpi = 300)

data_A82 <- data_all %>%
  dplyr::filter(dataset != "Original") %>%
  dplyr::mutate(
    dataset = stats::relevel(factor(dataset), ref = "USA"),
    group   = factor(paste0("G", cost, ints, dem)),
    group   = stats::relevel(group, ref = "G000")
  )

fit_A82 <- MASS::polr(
  strike5 ~ ally * group + male + age + ed4 + dataset,
  data = data_A82, Hess = TRUE
)

td_A82 <- tidy_polr_H(fit_A82) %>% dplyr::select(term, estimate, std.error, p.value, N)

ctx_codes <- c("G001","G010","G011","G100","G101","G110","G111")

rows_A82 <- tibble::tibble(
  row_order = 1:(1 + 7 + 7 + 3 + 5 + 1),
  Effect = c(
    "Alliance (vs no alliance) [in G000]",
    paste0(ctx_codes, " [baseline vs G000]"),
    paste0("Alliance × ", ctx_codes),
    "Male (vs female)",
    "Age (in decades)",
    "College+ (vs less than college)",
    "Brazil (vs USA)",
    "China (vs USA)",
    "India (vs USA)",
    "Russia (vs USA)",
    "UK (vs USA)",
    "N"
  ),
  term = c(
    "ally",
    paste0("group", ctx_codes),
    paste0("ally:group", ctx_codes),
    "male","age","ed4",
    "datasetBrazil","datasetChina","datasetIndia","datasetRussia","datasetUK",
    NA_character_
  )
)

extract_one <- function(td, term_candidates) {
  for (t in term_candidates) {
    hit <- td %>% dplyr::filter(term == t)
    if (nrow(hit) == 1) return(hit[1, c("estimate","std.error","p.value","N")])
  }
  tibble::tibble(estimate = NA_real_, std.error = NA_real_, p.value = NA_real_, N = NA_real_)
}

A82_rows <- lapply(seq_len(nrow(rows_A82)), function(i) {
  if (rows_A82$Effect[i] == "N") {
    return(tibble::tibble(
      row_order = rows_A82$row_order[i],
      Effect = rows_A82$Effect[i],
      estimate = NA_real_, std.error = NA_real_, p.value = NA_real_,
      N = stats::nobs(fit_A82)
    ))
  }
  term_i <- rows_A82$term[i]
  term_candidates <- if (grepl("^ally:group", term_i)) {
    c(term_i, paste0(sub("^ally:group", "group", term_i), ":ally"))
  } else {
    c(term_i)
  }
  out <- extract_one(td_A82, term_candidates)
  tibble::tibble(row_order = rows_A82$row_order[i], Effect = rows_A82$Effect[i]) %>% dplyr::bind_cols(out)
})

A8_2_Table <- dplyr::bind_rows(A82_rows) %>% dplyr::arrange(row_order) %>% dplyr::select(-row_order)
utils::write.csv(A8_2_Table, "A8_2_Table_context_pooled.csv", row.names = FALSE)

rows_A83 <- tibble::tibble(
  row_order = c(1, 2:8, 9:15, 16, 17, 18, 19),
  Effect = c(
    "Alliance (vs no alliance) [in G000]",
    paste0("Alliance × ", ctx_codes, " (", c("HC, LS, DEM","HC, HS, NDEM","HC, HS, DEM","LC, LS, NDEM","LC, LS, DEM","LC, HS, NDEM","LC, HS, DEM"), ")"),
    paste0(ctx_codes, " (", c("HC, LS, DEM","HC, HS, NDEM","HC, HS, DEM","LC, LS, NDEM","LC, LS, DEM","LC, HS, NDEM","LC, HS, DEM"), ") [baseline vs G000]"),
    "Male (vs female)",
    "Age (in decades)",
    "College+ (vs less than college)",
    "N"
  ),
  term_type = c(
    "ally",
    paste0("ally:group", ctx_codes),
    paste0("group", ctx_codes),
    "male","age","ed4",
    NA_character_
  )
)

country_cols_A83 <- list()

for (cty in countries) {
  df_cty <- data_all %>%
    dplyr::filter(dataset == cty) %>%
    dplyr::mutate(
      group = factor(paste0("G", cost, ints, dem)),
      group = stats::relevel(group, ref = "G000")
    )
  
  fit_cty <- MASS::polr(
    strike5 ~ ally * group + male + age + ed4,
    data = df_cty, Hess = TRUE
  )
  
  td <- tidy_polr_H(fit_cty) %>% dplyr::select(term, estimate, std.error, p.value)
  N_cty <- stats::nobs(fit_cty)
  
  get_row <- function(term_key) {
    if (is.na(term_key)) return(tibble::tibble(estimate=NA_real_, std.error=NA_real_, p.value=NA_real_))
    cand <- if (grepl("^ally:group", term_key)) {
      c(term_key, paste0(sub("^ally:group","group", term_key), ":ally"))
    } else {
      c(term_key)
    }
    for (t in cand) {
      hit <- td %>% dplyr::filter(term == t)
      if (nrow(hit) == 1) return(hit %>% dplyr::select(estimate, std.error, p.value))
    }
    tibble::tibble(estimate=NA_real_, std.error=NA_real_, p.value=NA_real_)
  }
  
  block <- lapply(seq_len(nrow(rows_A83)), function(i) {
    if (rows_A83$Effect[i] == "N") {
      return(tibble::tibble(
        row_order = rows_A83$row_order[i],
        !!paste0(cty, "_estimate") := NA_real_,
        !!paste0(cty, "_std.error") := NA_real_,
        !!paste0(cty, "_p.value") := NA_real_,
        !!paste0(cty, "_N") := as.numeric(N_cty)
      ))
    }
    vals <- get_row(rows_A83$term_type[i])
    tibble::tibble(
      row_order = rows_A83$row_order[i],
      !!paste0(cty, "_estimate") := vals$estimate,
      !!paste0(cty, "_std.error") := vals$std.error,
      !!paste0(cty, "_p.value") := vals$p.value,
      !!paste0(cty, "_N") := NA_real_
    )
  })
  
  country_cols_A83[[cty]] <- dplyr::bind_rows(block)
}

A8_3_Table <- rows_A83 %>% dplyr::select(row_order, Effect)
for (cty in countries) {
  A8_3_Table <- A8_3_Table %>% dplyr::left_join(country_cols_A83[[cty]], by = "row_order")
}
A8_3_Table <- A8_3_Table %>% dplyr::arrange(row_order) %>% dplyr::select(-row_order)
utils::write.csv(A8_3_Table, "A8_3_Table_context_by_country.csv", row.names = FALSE)

# ============================================================
# A9 — US waves coefficients + Party3 plot
# ============================================================

data_all_num <- data_all %>%
  dplyr::mutate(
    ally = as.numeric(haven::zap_labels(ally)),
    male = as.numeric(haven::zap_labels(male)),
    age  = as.numeric(haven::zap_labels(age)),
    ed4  = as.numeric(haven::zap_labels(ed4)),
    attention_pass = as.numeric(haven::zap_labels(attention_pass))
  )

data_USwaves_full <- data_all_num %>%
  dplyr::filter(dataset %in% c("Original", "USA")) %>%
  dplyr::mutate(
    US_wave = ifelse(dataset == "Original", "2017", "2025"),
    US_wave = factor(US_wave, levels = c("2025", "2017"))
  )

fit_USwaves_full <- MASS::polr(
  strike5 ~ ally * US_wave + male + age + ed4,
  data = data_USwaves_full, Hess = TRUE
)

extract_USwave_terms <- function(td) {
  td %>%
    dplyr::filter(term %in% c("US_wave2017", "ally:US_wave2017", "US_wave2017:ally",
                              "ally", "male", "age", "ed4")) %>%
    dplyr::mutate(
      Predictor = dplyr::case_when(
        term == "US_wave2017" ~ "2017 (vs 2025)",
        term %in% c("ally:US_wave2017", "US_wave2017:ally") ~ "Alliance × 2017",
        term == "ally" ~ "Alliance (vs no alliance) [in 2025]",
        term == "male" ~ "Male (vs female)",
        term == "age"  ~ "Age (in decades)",
        term == "ed4"  ~ "College+ (vs less than college)",
        TRUE ~ term
      )
    ) %>%
    dplyr::select(Predictor, estimate, std.error, p.value, N, term)
}

A9_Table_full <- extract_USwave_terms(tidy_polr_H(fit_USwaves_full))
utils::write.csv(A9_Table_full, "A9_Table_USwaves_full.csv", row.names = FALSE)

data_USwaves_attn <- data_all_num %>%
  dplyr::filter(attention_pass == 1) %>%
  dplyr::filter(dataset %in% c("Original", "USA")) %>%
  dplyr::mutate(
    US_wave = ifelse(dataset == "Original", "2017", "2025"),
    US_wave = factor(US_wave, levels = c("2025", "2017"))
  )

fit_USwaves_attn <- MASS::polr(
  strike5 ~ ally * US_wave + male + age + ed4,
  data = data_USwaves_attn, Hess = TRUE
)

A9_Table_ACpass <- extract_USwave_terms(tidy_polr_H(fit_USwaves_attn))
utils::write.csv(A9_Table_ACpass, "A9_Table_USwaves_ACpass.csv", row.names = FALSE)

dat_party <- haven::read_sav("merged_USwaves_for_FigA9.sav") %>%
  dplyr::mutate(
    dataset = as.character(dataset),
    ally = as.numeric(haven::zap_labels(ally)),
    party = as.numeric(haven::zap_labels(party)),
    strike5_support = as.numeric(haven::zap_labels(strike5_support))
  )

party3_summary <- dat_party %>%
  dplyr::mutate(
    party3 = dplyr::case_when(
      party %in% 1:3 ~ "Democrat",
      party == 4     ~ "Independent",
      party %in% 5:7 ~ "Republican",
      TRUE ~ NA_character_
    ),
    party3 = factor(party3, levels = c("Democrat","Independent","Republican")),
    ally_f = factor(ally, levels = c(0,1), labels = c("No Alliance","Alliance"))
  ) %>%
  dplyr::filter(!is.na(party3), !is.na(ally_f), !is.na(strike5_support)) %>%
  dplyr::group_by(dataset, party3, ally_f) %>%
  dplyr::summarise(
    N = dplyr::n(),
    mean_support = mean(strike5_support),
    sd_support = sd(strike5_support),
    se_support = sd_support / sqrt(N),
    ci_low  = mean_support + qnorm(0.025) * se_support,
    ci_high = mean_support + qnorm(0.975) * se_support,
    .groups = "drop"
  )

utils::write.csv(party3_summary, "A9_Fig_data_party3_by_ally_and_wave.csv", row.names = FALSE)

ally_palette <- c("No Alliance" = "#444444", "Alliance" = "#D55E00")
pd <- position_dodge(width = 0.55)

A9_Fig <- ggplot(party3_summary,
                 aes(x = mean_support, y = party3, color = ally_f)) +
  geom_errorbarh(aes(xmin = ci_low, xmax = ci_high),
                 position = pd, height = 0.12) +
  geom_point(position = pd, size = 2) +
  facet_wrap(~ dataset, nrow = 1) +
  scale_x_continuous(
    limits = c(1, 5),
    breaks = 1:5,
    labels = c("1\nOppose\nstrongly", "", "3\nNeither favor\nnor oppose", "", "5\nFavor\nstrongly")
  ) +
  scale_color_manual(values = ally_palette, breaks = c("No Alliance", "Alliance")) +
  labs(
    title = "Appendix 9 Figure: Mean support by party × alliance × wave",
    x = NULL, y = NULL, color = NULL
  ) +
  forest_style() +
  theme(
    legend.position = "bottom",
    panel.spacing.x = grid::unit(2.2, "lines"),
    axis.text.x = ggplot2::element_text(vjust = 1),
    plot.margin = grid::unit(c(5.5, 5.5, 10, 5.5), "pt")
  )

ggsave("A9_Fig_party3_by_ally_and_wave.png", A9_Fig, width = 10.5, height = 5.2, dpi = 300)

# ============================================================
# A10 — Analytical power curves (Wald)
# ============================================================

find_ally_interaction_term <- function(fit) {
  nm <- names(stats::coef(fit))
  it <- nm[nm != "ally" & grepl("ally", nm) & grepl(":", nm)]
  if (length(it) == 0) return(NA_character_)
  it[1]
}

get_polr_se <- function(fit, term_name) {
  cs <- coef(summary(fit))
  td <- as.data.frame(cs)
  td$term <- rownames(cs)
  rownames(td) <- NULL
  td <- td[!grepl("\\|", td$term), , drop = FALSE]
  td$`Std. Error`[match(term_name, td$term)]
}

power_curve_polr_analytic <- function(fit, term_name,
                                      effect_grid = seq(0, 0.6, by = 0.05),
                                      alpha = 0.05) {
  se <- get_polr_se(fit, term_name)
  zcrit <- stats::qnorm(1 - alpha/2)
  
  power_fn <- function(beta_true) {
    delta <- beta_true / se
    stats::pnorm(-zcrit - delta) + (1 - stats::pnorm(zcrit - delta))
  }
  
  pow <- vapply(effect_grid, power_fn, numeric(1))
  
  dplyr::tibble(
    beta  = effect_grid,
    OR    = exp(effect_grid),
    power = pow,
    term  = term_name,
    SE    = se,
    N     = stats::nobs(fit)
  )
}

term_H1 <- "ally"
term_H2 <- find_ally_interaction_term(fit_H2)
term_H3 <- find_ally_interaction_term(fit_H3)
term_H4 <- find_ally_interaction_term(fit_H4)

grid_beta <- seq(0, 0.6, by = 0.05)

pow_H1 <- power_curve_polr_analytic(fit_H1_pooled, term_H1, effect_grid = grid_beta) %>%
  dplyr::mutate(Hypothesis = "H1 (ally main effect)")

pow_H2 <- power_curve_polr_analytic(fit_H2, term_H2, effect_grid = grid_beta) %>%
  dplyr::mutate(Hypothesis = "H2 (ally×NATO interaction)")

pow_H3 <- power_curve_polr_analytic(fit_H3, term_H3, effect_grid = grid_beta) %>%
  dplyr::mutate(Hypothesis = "H3 (ally×Russia vs CIB interaction)")

pow_H4 <- power_curve_polr_analytic(fit_H4, term_H4, effect_grid = grid_beta) %>%
  dplyr::mutate(Hypothesis = "H4 (ally×USA vs UK interaction)")

A10_power_grid <- dplyr::bind_rows(pow_H1, pow_H2, pow_H3, pow_H4)
utils::write.csv(A10_power_grid, "A10_Table_power_curve_grid_ANALYTICAL.csv", row.names = FALSE)

A10_Fig <- ggplot2::ggplot(A10_power_grid, ggplot2::aes(x = beta, y = power, color = Hypothesis)) +
  ggplot2::geom_line(linewidth = 1) +
  ggplot2::geom_point(size = 1.5) +
  ggplot2::geom_hline(yintercept = 0.80, linetype = "dashed") +
  ggplot2::scale_x_continuous(breaks = seq(0, 0.6, by = 0.05), limits = c(0, 0.6)) +
  ggplot2::labs(
    title = "Appendix 10 Figure: Analytical power curves (Wald, p < .05)",
    x = "True effect size (beta, log-odds)",
    y = "Power",
    color = NULL
  ) +
  theme_pub +
  ggplot2::theme(legend.position = "bottom") +
  ggplot2::guides(color = ggplot2::guide_legend(nrow = 2, byrow = TRUE))

ggsave("A10_Fig_power_curves_ANALYTICAL.png", A10_Fig, width = 9, height = 5.5, dpi = 300)

A10_MDE80 <- A10_power_grid %>%
  dplyr::group_by(Hypothesis) %>%
  dplyr::summarise(
    SE       = dplyr::first(SE),
    N        = dplyr::first(N),
    MDE_beta = suppressWarnings(min(beta[power >= 0.80], na.rm = TRUE)),
    MDE_OR   = exp(MDE_beta),
    .groups  = "drop"
  )

utils::write.csv(A10_MDE80, "A10_Table_MDE_80pct_grid_ANALYTICAL.csv", row.names = FALSE)

# ============================================================
# A11 — Random-effects meta-analysis of ally effect (country models)
# ============================================================

data_pooled_meta <- data_all %>%
  dplyr::filter(dataset != "Original") %>%
  dplyr::mutate(
    dataset = factor(as.character(dataset)),
    strike5 = factor(as.character(strike5),
                     levels = as.character(0:4),
                     ordered = TRUE)
  )

tidy_polr_wald_meta <- function(fit) {
  broom::tidy(fit, conf.int = TRUE, conf.method = "Wald") %>%
    dplyr::filter(!grepl("\\|", term))
}

extract_ally_row_meta <- function(fit, label) {
  td <- tidy_polr_wald_meta(fit) %>% dplyr::filter(term == "ally")
  tibble::tibble(
    dataset   = label,
    estimate  = td$estimate,
    se        = td$std.error,
    conf.low  = td$conf.low,
    conf.high = td$conf.high,
    n         = stats::nobs(fit)
  )
}

meta_countries <- sort(unique(as.character(data_pooled_meta$dataset)))

ally_country <- lapply(meta_countries, function(cty) {
  df_cty <- data_pooled_meta %>% dplyr::filter(dataset == cty)
  
  fit_cty <- MASS::polr(
    strike5 ~ ally + male + age + ed4,
    data = df_cty, Hess = TRUE
  )
  
  extract_ally_row_meta(fit_cty, cty)
}) %>% dplyr::bind_rows()

utils::write.csv(ally_country, "A11_Table_meta_inputs_ally_country_estimates.csv", row.names = FALSE)

meta_re <- metafor::rma(
  yi  = estimate,
  sei = se,
  data = ally_country,
  method = "REML"
)

A11_meta_summary <- data.frame(
  model = "Random-effects (REML)",
  k = meta_re$k,
  overall_estimate = as.numeric(meta_re$b),
  overall_se = meta_re$se,
  overall_ci_low = meta_re$ci.lb,
  overall_ci_high = meta_re$ci.ub,
  tau2 = meta_re$tau2,
  I2 = meta_re$I2,
  H2 = meta_re$H2,
  Q = meta_re$QE,
  Q_df = meta_re$k - 1,
  Q_p = meta_re$QEp
)

utils::write.csv(A11_meta_summary, "A11_Table_meta_random_effects_summary.csv", row.names = FALSE)

A11_overall_or <- data.frame(
  OR = exp(as.numeric(meta_re$b)),
  OR_ci_low = exp(meta_re$ci.lb),
  OR_ci_high = exp(meta_re$ci.ub)
)
utils::write.csv(A11_overall_or, "A11_Table_meta_overall_OR.csv", row.names = FALSE)

meta_row <- tibble::tibble(
  dataset   = "Overall (RE meta)",
  estimate  = as.numeric(meta_re$b),
  se        = meta_re$se,
  conf.low  = meta_re$ci.lb,
  conf.high = meta_re$ci.ub,
  n         = NA_integer_
)

plot_df <- dplyr::bind_rows(meta_row, ally_country) %>%
  dplyr::mutate(dataset = factor(dataset, levels = c(meta_countries, "Overall (RE meta)")))

A11_Fig <- ggplot(plot_df, aes(y = dataset, x = estimate)) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), height = 0.15) +
  geom_point(aes(shape = dataset == "Overall (RE meta)"), size = 2.6) +
  scale_shape_manual(values = c(`TRUE` = 18, `FALSE` = 16), guide = "none") +
  labs(
    title = "Appendix 11 Figure: Random-effects meta-analysis of alliance effect",
    x = "Ally coefficient (log-odds, ordinal logit)",
    y = NULL
  ) +
  theme_pub

ggsave("A11_Fig_meta_ally_forestplot.png", A11_Fig, width = 7.8, height = 5.8, dpi = 300)


############################################################
## A7 — Mediation (product-of-coefficients)
############################################################

suppressPackageStartupMessages({
  library(emmeans)
  library(purrr)
})

emmeans::emm_options(disable = TRUE)

#set to 100 for testing the code, change to 1000 to replicate figures in the Appendix
B_boot <- 100

# ----------------------------
# Helpers
# ----------------------------

prep_mediation_df <- function(df) {
  df %>%
    dplyr::mutate(
      strike2  = as.numeric(haven::zap_labels(strike2)),
      r_hm     = as.numeric(haven::zap_labels(r_hm)),
      m_oblig  = as.numeric(haven::zap_labels(m_oblig)),
      ally     = as.numeric(haven::zap_labels(ally)),
      dem      = as.numeric(haven::zap_labels(dem)),
      ints     = as.numeric(haven::zap_labels(ints)),
      cost     = as.numeric(haven::zap_labels(cost)),
      male     = as.numeric(haven::zap_labels(male)),
      age      = as.numeric(haven::zap_labels(age)),
      ed4      = as.numeric(haven::zap_labels(ed4)),
      ally_f    = factor(ally),
      dem_f     = factor(dem),
      ints_f    = factor(ints),
      cost_f    = factor(cost),
      dataset_f = factor(dataset)
    ) %>%
    tidyr::drop_na(
      strike2, r_hm, m_oblig,
      ally_f, dem_f, ints_f, cost_f,
      male, age, ed4, dataset_f
    )
}

get_ally_effect <- function(model) {
  emm <- suppressMessages(
    emmeans::emmeans(model, specs = ~ ally_f, weights = "equal")
  )
  emmtab <- as.data.frame(emm)
  lv <- as.character(emmtab$ally_f)
  if (!all(c("0", "1") %in% lv)) return(NA_real_)
  emmtab$emmean[lv == "1"] - emmtab$emmean[lv == "0"]
}

med_product_once <- function(data, treat_term) {
  mediators <- c("r_hm", "m_oblig")
  controls  <- c("male", "age", "ed4")
  controls_term <- paste(controls, collapse = " + ")
  
  f_y <- as.formula(paste("strike2 ~", treat_term, "+", controls_term))
  f_y_m <- as.formula(paste("strike2 ~", paste(mediators, collapse = " + "), "+", treat_term, "+", controls_term))
  
  mod_y <- lm(f_y, data = data)
  tau   <- get_ally_effect(mod_y)
  
  a_paths <- purrr::map_dbl(mediators, function(m) {
    f_m <- as.formula(paste(m, "~", treat_term, "+", controls_term))
    mod_m <- lm(f_m, data = data)
    get_ally_effect(mod_m)
  })
  names(a_paths) <- mediators
  
  mod_y_m <- lm(f_y_m, data = data)
  b_paths <- stats::coef(mod_y_m)[mediators]
  
  ab <- a_paths * b_paths
  other_ab <- tau - sum(ab)
  
  shares <- 100 * ab / tau
  other_share <- 100 * other_ab / tau
  
  list(
    tau = as.numeric(tau),
    a   = a_paths,
    b   = b_paths,
    ab  = c(ab, Other = other_ab),
    shares = c(shares, Other = other_share)
  )
}

run_product_mediation_generic <- function(dfc, treat_term, B = 1000, seed = 54321) {
  N <- nrow(dfc)
  
  if (N < 30L) {
    paths <- tibble::tibble(
      Mediator      = c("r_hm", "m_oblig", "Other"),
      a             = NA_real_, a_low = NA_real_, a_high = NA_real_, a_pvalue = NA_real_,
      b             = NA_real_,
      ab            = NA_real_, ab_low = NA_real_, ab_high = NA_real_, ab_pvalue = NA_real_,
      Total_Effect  = NA_real_, Total_low = NA_real_, Total_high = NA_real_,
      Share_pct     = NA_real_, Share_low = NA_real_, Share_high = NA_real_,
      N             = N
    )
    diffs <- tibble::tibble(
      contrast     = c("ally_on_m_oblig_minus_r_hm", "indirect_m_oblig_minus_r_hm"),
      quantity     = c("a_path_effect", "ab_indirect_effect"),
      estimate     = NA_real_, low = NA_real_, high = NA_real_, pvalue = NA_real_,
      Total_Effect = NA_real_,
      N            = N
    )
    return(list(paths = paths, diffs = diffs))
  }
  
  point <- med_product_once(dfc, treat_term)
  share_point <- point$shares
  ab_point    <- point$ab
  a_point     <- point$a
  b_point     <- point$b
  tau_point   <- point$tau
  
  set.seed(seed)
  boot_list <- replicate(B, {
    idx <- sample.int(N, N, replace = TRUE)
    med_product_once(dfc[idx, , drop = FALSE], treat_term)
  }, simplify = FALSE)
  
  boot_share <- do.call(rbind, lapply(boot_list, `[[`, "shares"))
  boot_ab    <- do.call(rbind, lapply(boot_list, `[[`, "ab"))
  boot_a     <- do.call(rbind, lapply(boot_list, `[[`, "a"))
  boot_tau   <- vapply(boot_list, `[[`, numeric(1), "tau")
  
  colnames(boot_share) <- names(share_point)
  colnames(boot_ab)    <- names(ab_point)
  colnames(boot_a)     <- names(a_point)
  
  ci_share_low  <- apply(boot_share, 2, stats::quantile, probs = 0.025, na.rm = TRUE)
  ci_share_high <- apply(boot_share, 2, stats::quantile, probs = 0.975, na.rm = TRUE)
  ci_ab_low     <- apply(boot_ab,    2, stats::quantile, probs = 0.025, na.rm = TRUE)
  ci_ab_high    <- apply(boot_ab,    2, stats::quantile, probs = 0.975, na.rm = TRUE)
  ci_a_low      <- apply(boot_a,     2, stats::quantile, probs = 0.025, na.rm = TRUE)
  ci_a_high     <- apply(boot_a,     2, stats::quantile, probs = 0.975, na.rm = TRUE)
  ci_tau        <- stats::quantile(boot_tau, probs = c(0.025, 0.975), na.rm = TRUE)
  
  mediators_full <- names(share_point)  # r_hm, m_oblig, Other
  
  p_a  <- setNames(rep(NA_real_, length(mediators_full)), mediators_full)
  p_ab <- setNames(rep(NA_real_, length(mediators_full)), mediators_full)
  
  for (med in c("r_hm", "m_oblig")) {
    a_boot  <- boot_a[, med]
    ab_boot <- boot_ab[, med]
    p_a[med] <- 2 * min(mean(a_boot  >= 0, na.rm = TRUE), mean(a_boot  <= 0, na.rm = TRUE))
    p_ab[med] <- 2 * min(mean(ab_boot >= 0, na.rm = TRUE), mean(ab_boot <= 0, na.rm = TRUE))
  }
  
  paths <- tibble::tibble(
    Mediator      = mediators_full,
    a             = c(a_point, NA_real_),
    a_low         = c(ci_a_low, NA_real_),
    a_high        = c(ci_a_high, NA_real_),
    a_pvalue      = p_a[mediators_full],
    b             = c(b_point, NA_real_),
    ab            = ab_point[mediators_full],
    ab_low        = ci_ab_low[mediators_full],
    ab_high       = ci_ab_high[mediators_full],
    ab_pvalue     = p_ab[mediators_full],
    Total_Effect  = tau_point,
    Total_low     = ci_tau[1],
    Total_high    = ci_tau[2],
    Share_pct     = share_point[mediators_full],
    Share_low     = ci_share_low[mediators_full],
    Share_high    = ci_share_high[mediators_full],
    N             = N
  )
  
  diff_a_point  <- a_point["m_oblig"] - a_point["r_hm"]
  diff_ab_point <- ab_point["m_oblig"] - ab_point["r_hm"]
  
  diff_a_boot  <- boot_a[, "m_oblig"] - boot_a[, "r_hm"]
  diff_ab_boot <- boot_ab[, "m_oblig"] - boot_ab[, "r_hm"]
  
  diffs <- tibble::tibble(
    contrast     = c("ally_on_m_oblig_minus_r_hm", "indirect_m_oblig_minus_r_hm"),
    quantity     = c("a_path_effect", "ab_indirect_effect"),
    estimate     = c(diff_a_point, diff_ab_point),
    low          = c(stats::quantile(diff_a_boot,  0.025, na.rm = TRUE),
                     stats::quantile(diff_ab_boot, 0.025, na.rm = TRUE)),
    high         = c(stats::quantile(diff_a_boot,  0.975, na.rm = TRUE),
                     stats::quantile(diff_ab_boot, 0.975, na.rm = TRUE)),
    pvalue       = c(
      2 * min(mean(diff_a_boot  >= 0, na.rm = TRUE), mean(diff_a_boot  <= 0, na.rm = TRUE)),
      2 * min(mean(diff_ab_boot >= 0, na.rm = TRUE), mean(diff_ab_boot <= 0, na.rm = TRUE))
    ),
    Total_Effect = tau_point,
    N            = N
  )
  
  list(paths = paths, diffs = diffs)
}

run_product_mediation_country <- function(df_country, B = 1000, seed = 54321) {
  dfc <- prep_mediation_df(df_country)
  treat_term <- "ally_f * dem_f * ints_f * cost_f"
  run_product_mediation_generic(dfc, treat_term, B = B, seed = seed)
}

run_product_mediation_pooled <- function(df_all, B = 1000, seed = 54321) {
  dfc <- prep_mediation_df(df_all)
  treat_term <- "(ally_f + dem_f + ints_f + cost_f) * dataset_f"
  run_product_mediation_generic(dfc, treat_term, B = B, seed = seed)
}

# ----------------------------
# Run models
# ----------------------------

datasets_vec <- sort(unique(data_all$dataset))

med_paths_list <- list()
med_diffs_list <- list()

# Pooled (exclude Original)
res_pooled <- run_product_mediation_pooled(
  df_all = data_all %>% dplyr::filter(dataset != "Original"),
  B = B_boot, seed = 54321
)

A7_pooled_paths <- res_pooled$paths %>% dplyr::mutate(dataset = "Pooled", .before = 1)
A7_pooled_diffs <- res_pooled$diffs %>% dplyr::mutate(dataset = "Pooled", .before = 1)

utils::write.csv(A7_pooled_paths, "A7_Table_pooled_paths.csv", row.names = FALSE)
utils::write.csv(A7_pooled_diffs, "A7_Table_pooled_diffs.csv", row.names = FALSE)

med_paths_list[["Pooled"]] <- A7_pooled_paths
med_diffs_list[["Pooled"]] <- A7_pooled_diffs

# Per-dataset (includes Original)
for (ds in datasets_vec) {
  df_ds <- data_all %>% dplyr::filter(dataset == ds)
  res_ds <- run_product_mediation_country(df_ds, B = B_boot, seed = 54321)
  
  med_paths_list[[ds]] <- res_ds$paths %>% dplyr::mutate(dataset = ds, .before = 1)
  med_diffs_list[[ds]] <- res_ds$diffs %>% dplyr::mutate(dataset = ds, .before = 1)
}

A7_paths_all <- dplyr::bind_rows(med_paths_list)
A7_diffs_all <- dplyr::bind_rows(med_diffs_list)

utils::write.csv(A7_paths_all, "A7_Table_paths_all_datasets.csv", row.names = FALSE)
utils::write.csv(A7_diffs_all, "A7_Table_diffs_all_datasets.csv", row.names = FALSE)

# ----------------------------
# Figures
# ----------------------------

A7_pooled_plot <- A7_paths_all %>%
  dplyr::filter(dataset == "Pooled") %>%
  dplyr::mutate(
    Mediator_lab = dplyr::case_when(
      Mediator == "r_hm"    ~ "Reputation for Military Reliability",
      Mediator == "m_oblig" ~ "Moral Obligation",
      Mediator == "Other"   ~ "Other pathways",
      TRUE ~ Mediator
    ),
    Mediator_lab = factor(Mediator_lab, levels = c(
      "Other pathways",
      "Moral Obligation",
      "Reputation for Military Reliability"
    ))
  )

A7_Fig1 <- ggplot2::ggplot(A7_pooled_plot, aes(x = Share_pct, y = Mediator_lab)) +
  ggplot2::geom_vline(xintercept = 0, linetype = "dashed") +
  ggplot2::geom_errorbarh(aes(xmin = Share_low, xmax = Share_high), height = 0.15) +
  ggplot2::geom_point(size = 2) +
  ggplot2::geom_text(aes(label = round(Share_pct, 1)), vjust = -0.7, size = 3.3) +
  ggplot2::labs(
    x = "% of total alliance effect",
    y = NULL,
    title = "Appendix 7 Figure 2: Mechanism shares (pooled)"
  ) +
  theme_pub

ggplot2::ggsave("A7_Fig2_pooled_share.png", A7_Fig1, width = 6, height = 4, dpi = 300)

A7_pooled_a <- A7_paths_all %>%
  dplyr::filter(dataset == "Pooled", Mediator %in% c("r_hm", "m_oblig")) %>%
  dplyr::mutate(
    Mediator_lab = dplyr::case_when(
      Mediator == "r_hm"    ~ "Reputation for Military Reliability",
      Mediator == "m_oblig" ~ "Moral Obligation"
    ),
    Mediator_lab = factor(Mediator_lab, levels = c("Moral Obligation", "Reputation for Military Reliability"))
  )

A7_Fig2 <- ggplot2::ggplot(A7_pooled_a, aes(x = a, y = Mediator_lab)) +
  ggplot2::geom_vline(xintercept = 0, linetype = "dashed") +
  ggplot2::geom_errorbarh(aes(xmin = a_low, xmax = a_high), height = 0.15) +
  ggplot2::geom_point(size = 2) +
  ggplot2::geom_text(aes(label = round(a, 1)), vjust = -0.7, size = 3.3) +
  ggplot2::labs(
    x = "Effect of alliance on mediator (0–100 scale)",
    y = NULL,
    title = "Appendix 7 Figure 1: a-paths (pooled)"
  ) +
  theme_pub

ggplot2::ggsave("A7_Fig1_pooled_a.png", A7_Fig2, width = 6, height = 4, dpi = 300)

pos_j <- ggplot2::position_jitter(width = 0, height = 0.3, seed = 205)

A7_country_a <- A7_paths_all %>%
  dplyr::filter(dataset != "Pooled", dataset != "Original", Mediator %in% c("r_hm", "m_oblig")) %>%
  dplyr::mutate(
    Mediator_lab = dplyr::case_when(
      Mediator == "r_hm"    ~ "Reputation for Military Reliability",
      Mediator == "m_oblig" ~ "Moral Obligation"
    ),
    Mediator_lab = factor(Mediator_lab, levels = c("Moral Obligation", "Reputation for Military Reliability"))
  )

A7_Fig3 <- ggplot2::ggplot(A7_country_a, aes(x = a, y = Mediator_lab, color = dataset)) +
  ggplot2::geom_vline(xintercept = 0, linetype = "dashed") +
  ggplot2::geom_errorbarh(aes(xmin = a_low, xmax = a_high), position = pos_j, height = 0) +
  ggplot2::geom_point(position = pos_j, size = 2) +
  ggplot2::geom_text(aes(label = round(a, 1)), position = pos_j, vjust = -0.7, size = 2.8) +
  ggplot2::labs(
    x = "Effect of alliance on mediator (0–100 scale)",
    y = NULL,
    color = "Country",
    title = "Appendix 7 Figure 3: a-paths by country"
  ) +
  theme_pub +
  ggplot2::theme(legend.position = "bottom")

ggplot2::ggsave("A7_Fig3_country_a.png", A7_Fig3, width = 7, height = 4.5, dpi = 300)

A7_country_share <- A7_paths_all %>%
  dplyr::filter(dataset != "Pooled", dataset != "Original", Mediator %in% c("r_hm", "m_oblig")) %>%
  dplyr::mutate(
    Mediator_lab = dplyr::case_when(
      Mediator == "r_hm"    ~ "Reputation for Military Reliability",
      Mediator == "m_oblig" ~ "Moral Obligation"
    ),
    Mediator_lab = factor(Mediator_lab, levels = c("Moral Obligation", "Reputation for Military Reliability"))
  )

A7_Fig4 <- ggplot2::ggplot(A7_country_share, aes(x = Share_pct, y = Mediator_lab, color = dataset)) +
  ggplot2::geom_vline(xintercept = 0, linetype = "dashed") +
  ggplot2::geom_errorbarh(aes(xmin = Share_low, xmax = Share_high), position = pos_j, height = 0) +
  ggplot2::geom_point(position = pos_j, size = 2) +
  ggplot2::coord_cartesian(xlim = c(-50, 200)) +
  ggplot2::labs(
    x = "% of total alliance effect",
    y = NULL,
    color = "Country",
    title = "Appendix 7 Figure 4: mediator shares by country"
  ) +
  theme_pub +
  ggplot2::theme(legend.position = "bottom")

ggplot2::ggsave("A7_Fig4_country_share.png", A7_Fig4, width = 7, height = 4.5, dpi = 300)


############################################################
## MANUSCRIPT FIGURES 1–2 (STACKED)
## Uses objects created above:
##   - theme_pub, forest_style
##   - data_all, data_H1_pooled, fit_H1_pooled
##   - fit_H2, fit_H3, fit_H4
## Outputs:
##   - Fig1_H1_expected_then_coeffs_stacked.png
##   - Fig2_H2H4_expected_then_coeffs_stacked.png
############################################################

suppressPackageStartupMessages({
  library(patchwork)
})

# ----------------------------
# 0) Settings
# ----------------------------

R_pred <- 800

tick_labels_1to5_instr_fixed <- c(
  "Oppose strongly",              # 1
  "",                             # 2
  "Neither favor\nnor oppose",    # 3
  "",                             # 4
  "Favor strongly"                # 5
)

ally_palette <- c("No Alliance" = "#444444", "Alliance" = "#D55E00")

# ----------------------------
# 1) Prediction helpers (polr expected DV with ordered-threshold draws)
# ----------------------------

draw_theta_polr_ordered <- function(fit, R = 800, seed = 1, max_tries = 50) {
  set.seed(seed)
  
  b_hat <- stats::coef(fit)
  z_hat <- fit$zeta
  theta_hat <- c(b_hat, z_hat)
  
  V <- stats::vcov(fit)
  V <- V[names(theta_hat), names(theta_hat), drop = FALSE]
  
  keep <- NULL
  need <- R
  tries <- 0
  
  while (need > 0 && tries < max_tries) {
    tries <- tries + 1
    n_draw <- max(need * 2, 200)
    
    draws <- MASS::mvrnorm(n = n_draw, mu = theta_hat, Sigma = V)
    colnames(draws) <- names(theta_hat)
    
    z_draw <- draws[, names(z_hat), drop = FALSE]
    ok <- apply(z_draw, 1, function(z) all(diff(z) > 0))
    
    if (any(ok)) {
      good <- draws[ok, , drop = FALSE]
      take <- min(nrow(good), need)
      keep <- if (is.null(keep)) {
        good[seq_len(take), , drop = FALSE]
      } else {
        rbind(keep, good[seq_len(take), , drop = FALSE])
      }
      need <- R - nrow(keep)
    }
  }
  
  keep[seq_len(R), , drop = FALSE]
}

build_X_polr <- function(fit, newdata) {
  X_full <- stats::model.matrix(stats::delete.response(stats::terms(fit)), newdata)
  b_hat <- stats::coef(fit)
  X_full[, names(b_hat), drop = FALSE]
}

expected_0to4_draws <- function(X, beta_draw, zeta_draw) {
  eta <- X %*% t(beta_draw)  # N x R
  J <- ncol(zeta_draw)       # K-1
  
  cum_list <- vector("list", J)
  for (j in seq_len(J)) {
    cum_list[[j]] <- stats::plogis(sweep(-eta, 2, zeta_draw[, j], FUN = "+"))
  }
  
  probs <- vector("list", J + 1)
  probs[[1]] <- cum_list[[1]]
  if (J >= 2) {
    for (j in 2:J) probs[[j]] <- cum_list[[j]] - cum_list[[j - 1]]
  }
  probs[[J + 1]] <- 1 - cum_list[[J]]
  
  exp_mat <- 0 * probs[[1]]
  for (k in 2:(J + 1)) exp_mat <- exp_mat + (k - 1) * probs[[k]]
  
  colMeans(exp_mat)          # length R
}

expected_1to5_ci_for_scenario <- function(fit, newdata, theta_draws) {
  b_hat <- stats::coef(fit)
  z_hat <- fit$zeta
  
  b_draw <- theta_draws[, names(b_hat), drop = FALSE]
  z_draw <- theta_draws[, names(z_hat), drop = FALSE]
  
  X <- build_X_polr(fit, newdata)
  draws_0to4 <- expected_0to4_draws(X, b_draw, z_draw)
  ci <- stats::quantile(draws_0to4, probs = c(0.025, 0.975), na.rm = TRUE)
  
  dplyr::tibble(
    mean = mean(draws_0to4, na.rm = TRUE) + 1,
    low  = ci[[1]] + 1,
    high = ci[[2]] + 1
  )
}

expected_strike5_1to5_ci <- function(fit, R = 800, seed = 1) {
  df_cc <- stats::model.frame(fit)
  theta_draws <- draw_theta_polr_ordered(fit, R = R, seed = seed)
  
  d0 <- df_cc; d0$ally <- 0
  d1 <- df_cc; d1$ally <- 1
  
  out0 <- expected_1to5_ci_for_scenario(fit, d0, theta_draws)
  out1 <- expected_1to5_ci_for_scenario(fit, d1, theta_draws)
  
  dplyr::tibble(
    mean_no_alliance = out0$mean, low_no_alliance = out0$low, high_no_alliance = out0$high,
    mean_alliance    = out1$mean, low_alliance    = out1$low, high_alliance    = out1$high
  )
}

expected_by_group_and_ally_ci <- function(fit, group_var, group_levels, R = 800, seed = 1) {
  df_cc <- stats::model.frame(fit)
  theta_draws <- draw_theta_polr_ordered(fit, R = R, seed = seed)
  
  out_list <- lapply(group_levels, function(g) {
    d_base <- df_cc
    d_base[[group_var]] <- factor(g, levels = levels(df_cc[[group_var]]))
    
    d0 <- d_base; d0$ally <- 0
    d1 <- d_base; d1$ally <- 1
    
    o0 <- expected_1to5_ci_for_scenario(fit, d0, theta_draws) %>% dplyr::mutate(ally_status = "No Alliance")
    o1 <- expected_1to5_ci_for_scenario(fit, d1, theta_draws) %>% dplyr::mutate(ally_status = "Alliance")
    
    dplyr::bind_rows(o0, o1) %>% dplyr::mutate(subgroup = g)
  })
  
  dplyr::bind_rows(out_list)
}

# ----------------------------
# 2) Coefficient helpers (linear combination for polr slopes)
# ----------------------------

lincom_polr <- function(fit, L_named) {
  b <- stats::coef(fit)
  V <- stats::vcov(fit)
  
  slope_names <- names(b)
  Vb <- V[slope_names, slope_names, drop = FALSE]
  
  L <- rep(0, length(slope_names))
  names(L) <- slope_names
  common <- intersect(names(L_named), slope_names)
  L[common] <- L_named[common]
  
  est <- sum(L * b)
  se  <- sqrt(as.numeric(t(L) %*% Vb %*% L))
  ci  <- est + c(-1, 1) * stats::qnorm(0.975) * se
  
  dplyr::tibble(
    estimate = est,
    conf.low = ci[[1]],
    conf.high = ci[[2]],
    se = se
  )
}

find_interaction_term <- function(fit, group_var) {
  nm <- names(stats::coef(fit))
  it <- nm[grepl(paste0("^ally:", group_var), nm) | (grepl(paste0("^", group_var), nm) & grepl(":ally$", nm))]
  if (length(it) == 0) return(NA_character_)
  it[1]
}

extract_first_ally_interaction <- function(fit, Hypothesis) {
  td <- broom::tidy(fit, conf.int = TRUE, conf.method = "Wald") %>%
    dplyr::filter(!grepl("\\|", term)) %>%
    dplyr::filter(grepl("^ally:", term) | grepl(":ally$", term)) %>%
    dplyr::slice(1)
  
  dplyr::tibble(
    Hypothesis  = Hypothesis,
    Subsample   = "Interaction",
    estimate    = td$estimate,
    conf.low    = td$conf.low,
    conf.high   = td$conf.high,
    effect_type = "Interaction"
  )
}

# ============================================================
# FIGURE 1 (H1): Expected DV stacked over ally coefficients
# ============================================================

# Country-specific H1 fits (same spec as Appendix A5.2)
fits_H1_country <- setNames(vector("list", length(countries)), countries)
for (cty in countries) {
  df_cty <- data_H1_pooled %>% dplyr::filter(as.character(dataset) == cty)
  fits_H1_country[[cty]] <- MASS::polr(strike5 ~ ally + male + age + ed4, data = df_cty, Hess = TRUE)
}

country_order_topdown <- c("Pooled","USA","UK","Russia","India","China","Brazil")

# Expected DV (pooled + countries)
pred_pooled <- expected_strike5_1to5_ci(fit_H1_pooled, R = R_pred, seed = 150) %>%
  dplyr::mutate(dataset = "Pooled", .before = 1)

pred_countries <- dplyr::bind_rows(lapply(names(fits_H1_country), function(cty) {
  expected_strike5_1to5_ci(fits_H1_country[[cty]], R = R_pred, seed = 100 + match(cty, names(fits_H1_country))) %>%
    dplyr::mutate(dataset = cty, .before = 1)
}))

pred_H1 <- dplyr::bind_rows(pred_pooled, pred_countries) %>%
  dplyr::mutate(dataset = factor(dataset, levels = rev(country_order_topdown)))

pred_H1_long <- pred_H1 %>%
  tidyr::pivot_longer(
    cols = c(mean_no_alliance, low_no_alliance, high_no_alliance,
             mean_alliance,    low_alliance,    high_alliance),
    names_to = c(".value", "ally_status"),
    names_pattern = "(mean|low|high)_(no_alliance|alliance)"
  ) %>%
  dplyr::mutate(
    ally_status = factor(ally_status, levels = c("no_alliance", "alliance"),
                         labels = c("No Alliance", "Alliance")),
    dataset = factor(dataset, levels = levels(pred_H1$dataset))
  )

pdodge_h1 <- ggplot2::position_dodge(width = 0.55)

p_H1_pred_forest <- ggplot2::ggplot(
  pred_H1_long,
  ggplot2::aes(y = dataset, x = mean, color = ally_status)
) +
  ggplot2::geom_errorbarh(
    ggplot2::aes(xmin = low, xmax = high),
    position = pdodge_h1,
    height = 0.12
  ) +
  ggplot2::geom_point(position = pdodge_h1, size = 2) +
  ggplot2::scale_x_continuous(
    limits = c(1, 5),
    breaks = 1:5,
    labels = tick_labels_1to5_instr_fixed
  ) +
  ggplot2::scale_color_manual(values = ally_palette, breaks = c("No Alliance", "Alliance")) +
  ggplot2::labs(
    title = NULL,
    x = "Expected support for strike",
    y = NULL,
    color = NULL
  ) +
  forest_style() +
  ggplot2::theme(legend.position = "bottom")

# Ally coefficients (pooled + countries)
ally_rows <- dplyr::bind_rows(
  broom::tidy(fit_H1_pooled, conf.int = TRUE, conf.method = "Wald") %>%
    dplyr::filter(!grepl("\\|", term), term == "ally") %>%
    dplyr::transmute(dataset = "Pooled", estimate, conf.low, conf.high),
  dplyr::bind_rows(lapply(names(fits_H1_country), function(cty) {
    broom::tidy(fits_H1_country[[cty]], conf.int = TRUE, conf.method = "Wald") %>%
      dplyr::filter(!grepl("\\|", term), term == "ally") %>%
      dplyr::transmute(dataset = cty, estimate, conf.low, conf.high)
  }))
) %>%
  dplyr::mutate(dataset = factor(dataset, levels = rev(country_order_topdown)))

p_ally_full <- ggplot2::ggplot(ally_rows, ggplot2::aes(y = dataset, x = estimate)) +
  ggplot2::geom_vline(xintercept = 0, linetype = "dashed") +
  ggplot2::geom_errorbarh(ggplot2::aes(xmin = conf.low, xmax = conf.high), height = 0.12) +
  ggplot2::geom_point(size = 2) +
  ggplot2::labs(
    title = NULL,
    x = "Alliance effect (ordinal logit, log-odds)",
    y = NULL
  ) +
  forest_style()

Fig1_H1_stacked <- p_H1_pred_forest / p_ally_full +
  patchwork::plot_layout(ncol = 1, heights = c(1, 1))

ggplot2::ggsave(
  "Fig1_H1_expected_then_coeffs_stacked.png",
  Fig1_H1_stacked,
  width = 9.5, height = 11.5, dpi = 300
)

# ============================================================
# FIGURE 2 (H2–H4): Expected DV stacked over implied coefficients
# ============================================================

# Expected DV (by subgroup × alliance)
pred_H2_long <- expected_by_group_and_ally_ci(
  fit_H2, group_var = "NATO",
  group_levels = levels(stats::model.frame(fit_H2)$NATO),
  R = R_pred, seed = 201
) %>% dplyr::mutate(Hypothesis = "H2")

pred_H3_long <- expected_by_group_and_ally_ci(
  fit_H3, group_var = "Rus_vs_CIB",
  group_levels = levels(stats::model.frame(fit_H3)$Rus_vs_CIB),
  R = R_pred, seed = 301
) %>%
  dplyr::mutate(
    Hypothesis = "H3",
    subgroup = ifelse(subgroup == "China/India/Brazil", "China / India / Brazil", subgroup)
  )

pred_H4_long <- expected_by_group_and_ally_ci(
  fit_H4, group_var = "US_vs_UK",
  group_levels = levels(stats::model.frame(fit_H4)$US_vs_UK),
  R = R_pred, seed = 401
) %>% dplyr::mutate(Hypothesis = "H4")

pred_H2H4_long <- dplyr::bind_rows(pred_H2_long, pred_H3_long, pred_H4_long) %>%
  dplyr::mutate(
    ally_status = factor(ally_status, levels = c("No Alliance", "Alliance")),
    Hypothesis  = factor(Hypothesis, levels = c("H2", "H3", "H4")),
    subgroup    = factor(as.character(subgroup),
                         levels = c("Non-NATO","NATO",
                                    "China / India / Brazil","Russia",
                                    "USA","UK"))
  )

# Implied ally coefficients within subgroups (from interaction models)
it_H2 <- find_interaction_term(fit_H2, "NATO")
it_H3 <- find_interaction_term(fit_H3, "Rus_vs_CIB")
it_H4 <- find_interaction_term(fit_H4, "US_vs_UK")

coef_H2 <- dplyr::bind_rows(
  lincom_polr(fit_H2, c(ally = 1)) %>% dplyr::mutate(Subsample = "Non-NATO", Hypothesis = "H2"),
  lincom_polr(fit_H2, stats::setNames(c(1, 1), c("ally", it_H2))) %>% dplyr::mutate(Subsample = "NATO", Hypothesis = "H2")
)

coef_H3 <- dplyr::bind_rows(
  lincom_polr(fit_H3, c(ally = 1)) %>% dplyr::mutate(Subsample = "China / India / Brazil", Hypothesis = "H3"),
  lincom_polr(fit_H3, stats::setNames(c(1, 1), c("ally", it_H3))) %>% dplyr::mutate(Subsample = "Russia", Hypothesis = "H3")
)

coef_H4 <- dplyr::bind_rows(
  lincom_polr(fit_H4, c(ally = 1)) %>% dplyr::mutate(Subsample = "USA", Hypothesis = "H4"),
  lincom_polr(fit_H4, stats::setNames(c(1, 1), c("ally", it_H4))) %>% dplyr::mutate(Subsample = "UK", Hypothesis = "H4")
)

coef_subsamples <- dplyr::bind_rows(coef_H2, coef_H3, coef_H4) %>%
  dplyr::mutate(
    Hypothesis = factor(Hypothesis, levels = c("H2","H3","H4")),
    Subsample  = factor(Subsample, levels = c("Non-NATO","NATO",
                                              "China / India / Brazil","Russia",
                                              "USA","UK"))
  )

# Add interaction coefficient rows (as their own middle row per hypothesis)
int_rows <- dplyr::bind_rows(
  extract_first_ally_interaction(fit_H2, "H2"),
  extract_first_ally_interaction(fit_H3, "H3"),
  extract_first_ally_interaction(fit_H4, "H4")
)

coef_rows <- coef_subsamples %>%
  dplyr::transmute(
    Hypothesis  = as.character(Hypothesis),
    Subsample   = as.character(Subsample),
    estimate,
    conf.low,
    conf.high,
    effect_type = "Alliance effect"
  )

coef_with_int <- dplyr::bind_rows(coef_rows, int_rows) %>%
  dplyr::mutate(
    Hypothesis  = factor(Hypothesis, levels = c("H2", "H3", "H4")),
    Subsample   = as.character(Subsample),
    y_key       = paste0(as.character(Hypothesis), "__", Subsample),
    effect_type = factor(effect_type, levels = c("Alliance effect", "Interaction"))
  )

y_levels <- c(
  "H2__Non-NATO", "H2__Interaction", "H2__NATO",
  "H3__China / India / Brazil", "H3__Interaction", "H3__Russia",
  "H4__USA", "H4__Interaction", "H4__UK"
)

coef_with_int <- coef_with_int %>%
  dplyr::mutate(y_key = factor(y_key, levels = y_levels))

pdodge_h234 <- ggplot2::position_dodge(width = 0.55)

p_H2H4_exp_stacked <- ggplot2::ggplot(
  pred_H2H4_long,
  ggplot2::aes(y = subgroup, x = mean, color = ally_status)
) +
  ggplot2::geom_vline(xintercept = 1:5, color = "grey80", linewidth = 0.5) +
  ggplot2::geom_errorbarh(
    ggplot2::aes(xmin = low, xmax = high),
    position = pdodge_h234,
    height = 0.12
  ) +
  ggplot2::geom_point(position = pdodge_h234, size = 2) +
  ggplot2::facet_grid(
    Hypothesis ~ .,
    scales = "free_y",
    space  = "free_y",
    switch = "y"
  ) +
  ggplot2::scale_x_continuous(
    limits = c(1, 5),
    breaks = 1:5,
    labels = tick_labels_1to5_instr_fixed
  ) +
  ggplot2::scale_color_manual(values = ally_palette, breaks = c("No Alliance", "Alliance")) +
  ggplot2::labs(
    title = NULL,
    x = "Expected support for strike",
    y = NULL,
    color = NULL
  ) +
  forest_style() +
  ggplot2::theme(
    strip.placement   = "outside",
    strip.background  = ggplot2::element_blank(),
    strip.text.y.left = ggplot2::element_text(angle = 0, hjust = 0),
    panel.spacing.y   = grid::unit(2.2, "lines"),
    legend.position   = "bottom"
  )

p_H2H4_coef_stacked <- ggplot2::ggplot(
  coef_with_int,
  ggplot2::aes(y = y_key, x = estimate, color = effect_type)
) +
  ggplot2::geom_vline(xintercept = 0, linetype = "dashed") +
  ggplot2::geom_errorbarh(
    ggplot2::aes(xmin = conf.low, xmax = conf.high),
    height = 0.12
  ) +
  ggplot2::geom_point(size = 2) +
  ggplot2::facet_grid(
    Hypothesis ~ .,
    scales = "free_y",
    space  = "free_y",
    switch = "y"
  ) +
  ggplot2::scale_color_manual(
    values = c("Alliance effect" = "black", "Interaction" = "#2C7FB8"),
    breaks = c("Alliance effect", "Interaction")
  ) +
  ggplot2::scale_y_discrete(
    labels = function(x) {
      lab <- sub("^[^_]+__", "", x)
      ifelse(lab == "Interaction", "", lab)
    }
  ) +
  ggplot2::labs(
    title = NULL,
    x = "Alliance effect (ordinal logit, log-odds)",
    y = NULL,
    color = NULL
  ) +
  forest_style() +
  ggplot2::theme(
    strip.placement   = "outside",
    strip.background  = ggplot2::element_blank(),
    strip.text.y.left = ggplot2::element_text(angle = 0, hjust = 0),
    panel.spacing.y   = grid::unit(2.2, "lines"),
    legend.position   = "bottom"
  )

Fig2_H2H4_stacked <- p_H2H4_exp_stacked / p_H2H4_coef_stacked +
  patchwork::plot_layout(ncol = 1, heights = c(1, 1))

ggplot2::ggsave(
  "Fig2_H2H4_expected_then_coeffs_stacked.png",
  Fig2_H2H4_stacked,
  width = 9.5, height = 11.5, dpi = 300
)